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Summary 

The  solution  to  sets  of  linear  equations  is  fundamental  to  mathematical 
modeling.  The  choice  of  a  programming  language  for  the  algorithm  used  can 
critically  influence  the  efficiency  of  the  methodology.  The  viability  of  the 
MUMPS  programming  language  is  demonstrated  in  developing  an  extremely  conser¬ 
vative  and  accurate  routine  to  solve  linear  systems  using  a  complete  pivotal 
Gaussian  forward  elimination  strategy. 

The  technique  of  using  a  simple  hashing  scheme  to  avoid  post-multiplica¬ 
tive  matrices  to  identify  variable  domains  permits  tight,  concise  coding  while 
guaranteeing  maximum  accuracy,  given  the  conditions  of  the  problem.  Linear 
algebra  notation  combined  with  using  VAFM  in  a  relational  context  allows  the 
analysis  maximum  freedom  in  data  manipulation. 

3 

The  Gauss  algorithm  presented  in  this  paper  solves  any  n  by  n  system  in  n 

2 

steps.  Additional  solutions  demand  only  n  operations.  Complete  pivotal 
Gaussian  elimination  is  a  "perfect"  algorithm  in  a  mathematical  sense. 


OOPY 

.INSPECTED 


Accesion  For 

NTIS  CRA&I 
DTlC  TAB 
Unannounced 
Justification 


By . . 

Distribution/ 


Availability  Codes 


Avail  and/or 
Speeia! 


H 


COMPLETE  PIVOTAL  GAUSSIAN  ELIMINATION 
USING  THE  VETERANS'  ADMINISTRATION  FILE  MANAGER 


Dallas  R.  Hodgins 
Naval  Health  Research  Center 

Introduction 

The  installation  of  the  Naval  Occupational  Health  Information  Management 
System  (NOHIMS)  at  Naval  industrial  facilities  makes  possible  a  valuable  data 
resource  for  industrial  hygienists,  medical  researchers,  and  epidemiologists. 
The  Veterans'  Administration  File  Manager  (FM)  component  of  NOHIMS  manages 
this  data  linking  the  worker  with  his  environment.  Prior  to  implementation,  a 
need  to  augment  the  statistical  tools  of  FM  was  recognized.  The  first  phase 
of  this  effort  consisted  of  providing  researchers  with  basic  descriptive 
statistical  algorithms.  These  programs  were  written  in  MUMPS  (Massachusetts 
General  Hospital  Utility  Multi-Programming  System)  using  linear  algebra 
notation  and  embedded  in  FM  (also  written  in  MUMPS).  A  conscious  choice  was 
made  to  design  routines  to  efficiently  handle  small  data  samples,  realistic¬ 
ally  recognizing  the  magnitude  of  data  handled  and  the  response  time  require¬ 
ments  of  shipyard  industrial  hygienists.  Three  fundamental  issues,  however, 
remained  unresolved: 

(1)  Is  MUMPS  a  reasonable  choice  for  numerical  analysis  in  statistical 
work? 

(2)  Is  FM  the  right  choice  foi  statistical  data  structuring  and 
manipulation? 

(3)  Assuming  one  uses  MUMPS  and  FM,  what  is  the  best  logical  data 
structure? 

These  questions  are  addressed  in  this  paper  in  a  pragmatic  fashion  by 
developing  an  algorithm  to  solve  linear  equations.  Of  the  two  problems  cen¬ 
tral  to  abstract  linear  algebra,  solutions  to  Ax=b  and  Ax=\x,  the  first  has  a 
perfect  algorithm  (in  the  mathematical  sense)--Gaussian  elimination.  Indeed, 
the  basic  algorithm  of  numerical  lineal  algebra  is  Gaussian  elimination  with 
partial  pivoting.  This  paper  presents  a  complete  pivotal  strategy  using  a 
hashing  technique  to  eliminate  the  need  of  postmultiplication  by  a  permutation 
matrix  to  ascertain  the  original  domains  of  the  unknowns. 


The  overwhelming  power,  both  memory  and  processing,  of  current  machines 
has  led  to  the  use  of  "black  box"  algorithms  that  may  or  may  not  be  suitable 
to  the  problem  at  hand.  Numerical  techniques  designed  to  handle  very  large 
models  cannot  afford  the  luxury  of  the  conservativeness  of  a  complete  pivotal 
strategy  in  the  standard  decomposition  of  the  matrix  A.  In  keeping  with 
"Small  is  Beautiful",  the  numerical  model  presented  is  designed  for  problems 
encountered  by  medical  researchers,  industrial  hygienists,  and  statisticians 
that  are  dealing  with  a  handful  of  variables  and  at  most  a  few  hundred  cases 
— not  thousands.  These  people  want  to  know  if  their  mathematical  model  is 
stable  and  they  want  the  best  numerical  approach  to  the  solution.  They  are 
willing  to  sacrifice  machine  time  for  soundness  and  accuracy  in  their  re¬ 
search.  The  use  of  linear  algebra  and  the  "view"  of  the  data  afforded  by  FM, 
used  in  a  relational  sense,  coalesce  to  provide  a  powerful  and  practical 
environment  for  numerical  statistical  analysis. 

Materials  and  Methods 

This  work  was  done  using  Intersystems  MUMPS  M/VX  3.1  on  a  Digital  Equip¬ 
ment  Corporation  VAX  11/785  with  FM  version  17.07. 

The  mathematical  formulation  most  commonly  used  in  statistics  to  express 
relationships  between  two  or  more  variables  is  the  linear  equation  of  the  form 

Y'  =  X1  +  X2V1  +  X3V2  "•  XkVk  (1) 

where  Y'  is  to  be  predicted,  V^Vj,  ...  are  the  known  values  from  the 

observed  data  on  which  the  predictions  are  tn  be  based,  and  X^,X9  ...  X^  are 

the  numerical  coefficients  to  be  determined.  We  find  these  by  minimizing  the 

2 

sum  of  the  squares  (Y-Y' )  where  the  Y's  are  the  observed  values  and  the  Y 

primes  are  the  calculated  values.  This  is  done  by  taking  the  partial  deriva- 

2 

tives  of  the  expression  for  (Y-Y')  with  respect  to  as  many  unknown  coeffi¬ 
cients  Xj,  X2,  etc.  as  we  have  in  the  relationship. 

The  "normal"  equations  for  two  independent  variables  are: 


ZNX1 

+ 

(ZV^ 

+ 

(Zv2)x3 

ZY 

(ZvpXj 

+ 

(ZV1z)x2 

+ 

(Zv^^ 

ZVjY  (2) 

(Zv2)Xi 

+ 

(ZV1V2)X2 

+ 

(ZV22)X3 

-  2V 

where  N  is  the  number 

of 

observations 

1  » 

ZV1?  IV2  are 

the  sums  of  the  indepen 

dent  variables,  and  ZV 

1V2 

is  the  sum 

of 

the  product 

of  the  independent  varia 

bles  and  ZY  is  the  sum  of  the  inhomogeneous  terms. 
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More  concisely,  we  can  represent  this  system  by  Ax=b  where  A  is  the  square 
matrix  of  the  independent  variable  weights  of  the  overdetermined  system  in 
(1),  X  is  a  column  vector  representing  the  unknown  values  we  are  seeking,  and 
b  is  the  column  vector  of  inhomogeneous  terms.  Ax=b  has  a  solution  if  and 
only  if  b  lies  in  the  vector  space  spanned  by  the  column  vectors  of  A  (which 
we  made  orthogonal  to  the  error  vector  (y-y')  by  setting  the  partial  deriva¬ 
tions  to  zero,  hence  the  term  "normal").  The  solution  xsA'^b  will  be  demon¬ 
strated  using  a  complete  pivotal  forward  Gaussian  elimination  strategy.  The 
matrix  A  is  factored  into  lower  and  upper  triangular  matrices  L  and  U  such 
that  LU=A,  allowing  us  to  solve  for  A-^b,  not  A~^  itself.  The  original 
columns  of  the  unknowns  are  ascertained  at  completion  of  the  computation  by  a 
simple  hashing  technique  based  on  their  original  position  in  the  matrix.  This 
avoids  complicated  bookkeeping  needed  to  effect  the  post-multiplication  by 
permutation  matrices — the  price  of  which,  in  the  traditional  algebraic 
solution,  is  precisely  why  one  seldom,  if  ever,  finds  an  algorithm  that 
incorporates  complete  pivoting.  We  can  afford  the  conservatism  of  complete 
pivoting,  however,  as  we  are  expressly  manipulating  small  systems. 

Purely  hypothetically,  suppose  it  was  desired  to  predict  the  effect  on  the 
average  price  of  a  mobile  surgical  tent  of  the  number  of  operating  tables  and 
light  fixtures  to  illuminate  the  tables.  Using  FM  we  would  search  (using  the 
option  SEARCH  FILE  ENTRIES)  the  file  containing  fiscal  data  for  information 
relating  to  surgical  tents.  This  data  would  be  placed  in  the  global  ‘DIBT  for 
us  by  FM.  Using  the  TRANSFER  ENTRIES  option  we  transfer  the  following  data  to 
the  file  SURGTENT  (created  by  MODIFY  FILE  ATTRIBUTES  prior  to  the  transfer) 
from  ‘DIBT: 


Table  1.  Relation  SURGTENT: 


Number  of  Tables  (V^) 
3-Tuple  1  3 

2 

4 
2 

3 
2 

5 

4 


Number  of  Lights  (Vj) 
2 
1 
3 
1 
2 
2 
3 


Price  (Y),  Dollars 
33800 
29300 
38800 
29200 
34700 
29900 
43400 
37900 


3-Tuple  8 


Table  1  in  our  notation  is  X1(J,K)  where  XI  is  the  relation  SURGTENT,  J  is 
the  record  number  (the  primary  key — tuple  1,  tuple  2,  etc.)  and  K  is  the 
domain:  K=l=Number  of  Tables,  K=2=Number  of  Lights,  K=3=Price.  Therefore, 

Xl(7,3)  represents  the  table  (relation)  surgical  tent,  3-tuple  7  (set),  domain 
3  (price) — more  precisely — 43,400  dollars. 

Carefully  multiplying  and  summing,  we  obtain  the  last  line  in  Table  2: 

Table  2.  Independent  Variables 


V1 

V2 

Y 

V1Y 

¥  V1  V1V2 

V  z 

2 

3 

2 

33,800 

101,400 

67,600  9  6 

4 

2 

1 

29,300 

58,600 

29,300  4  2 

1 

4 

3 

38,000 

155,200 

116,400  16  12 

9 

2 

1 

29,200 

58,400 

29,200  4  2 

1 

3 

2 

34,700 

104,100 

69,400  9  6 

4 

2 

2 

29,900 

59,800 

59,800  4  4 

4 

5 

3 

43,400 

217,000 

130,200  25  15 

9 

4 _ 

2 

37,900 

151,600 

75,800  16  8 

4 

ZV1=25 

SV-16 

2)Y=277,000 

2^=906,100  2V2Y=577,700  2^=87  2^=55 

ZV22=36 

Inserting  the 

sums  into 

i  our  "normal"  equations  we  arrive  at 

8X  +  25X„ 

+  16X3  =  277000 

25XJ  +  87X2 

+  55X3  =  906100 

(3) 

16X.  +  55X„ 

+  36X3  =  577700 

the  "normal 

"  equations  for 

the  overdetermined  system  X1(J,K).  For  reasons  of 

exposi tion, 

a  file  NORMEQ 

(file  4003) 

was  created  to  hold  these 

values.  Part 

of  the  data 

dictionary  for 

NORMEQ  look 

s  like  this: 

CURRENT 

ENTRIES 

NAME:  EQ1 

XI:  8 

X2:  25 

X3:  16  X4:  277000 

NAME:  EQ2 

XI:  25 

X2:  87 

X3:  55  X4:  906100 

(4) 

NAME:  EQ3 

XI:  16 

X2:  55 

X3 :  36  X4:  577700 

which  does 

not  require  too 

much  imagination  to  visualize  as  the 

equations  in 

Let  us  digress  a  moment  to  familiarize  ourselves  with  the  terms  and  con¬ 
cepts  of  Gaussian  elimination  to  motivate  the  subsequent  discussion  of  the 

algorithm  for  the  multiple  regression  solution.  Consider  this  simple  system 

2 

of  three  equations  in  three  unknowns  (from  Strang  ): 


+  X_  +  2X. 


0  +  4X„ 


2XX  +  > 

which  can  be  expressed  as  Ax=b  where 


2X3  = 


and  b  = 


Simple  Elimination 

Displaying  our  problem  in  the  following  manner 


1  0  4 


=  -2 


we  proceed  as  follows: 

Step  1:  (a)  Divide  the  first  element  of  the  first  row  (the  first  pivot) 
into  the  first  element  of  the  second  row;  multiply  the  first  row  by  this  quo¬ 
tient  (equals  1)  and  subtract  the  first  row  from  the  second;  (b)  divide  the 
first  element  of  the  first  row  into  the  first  element  of  the  third  row;  multi¬ 
ply  the  first  row  by  this  quotient  (equals  2)  and  subtract  the  result  from  the 
third  row  obtaining: 


Step  2:  Divide  the  second  element  of  the  third  row  by  the  second  element 
of  the  second  row  (the  second  pivot);  multiply  the  second  row  by  this  quotient 
(equals  1)  and  subtract  the  second  row  from  the  third,  obtaining: 

1  1  2 j  fx  1  [l 

I  1 


0-1  2 
0  0  -8 


=  -3 


L  l  3  J  l  J 

This  completes  the  forward  elimination  allowing  us  to  solve  for  X3 
diately  (-8X^  =8,  X^  =  -1). 

Step  3:  Substituting  the  value  of  X^  back  into  row  two  we  obtain 

-X2  4  2(1)  =  -3  or  X2  =  1 

Back  substituting  the  values  of  X-^  and  X7  into  row  one  we  have 

X.  4  1  +  2(1)  -  1  or  X.  =  2. 


VVoV* 


Elimination  with  Partial  Pivoting 

Solving  the  same  system  with  a  partial  pivoting  strategy  we  simply  start 
the  forward  elimination  by  choosing  as  the  pivot  row  that  row  which  has  the 
largest  element  (in  an  absolute  sense)  in  the  pivot  column  and  repeat  the 


process  for  choosing  the  pivots  at  each  stage  of  the  elimination. 


where 


i  i 

2 

2 

1 

-2 

2 

1 

-2 

2 

1 

2 

A  = 

1  0 

4 

— ** 

1 

1 

2 

- 

0 

.5 

3 

0 

.5 

3 

2  1 

-2 

1 

0 

4 

0 

-.5 

5 

0 

0 

8 

i-2ir 

7 

r* 

i 

2 

0.5  3 

0  0  8 


2 

X- 
L  3 


-2.5 

-8 


again  yields  X 


or  X, 


2,  X2=l,  X3=-l 


Note  that  .5  and  -.5  are  of  the  same  absolute  magnitude  in  the  pivot  column  at 
the  second  stage  so  no  switching  was  necessary.  (The  operations  on  the  column 
vector  of  inhomogeneous  terms  are  not  shown  as  a  challenge  to  the  reader.) 
This  is  the  basic  numerical  algorithm  of  linear  algebra-- forward  elimination 
with  partial  pivoting. 


Elimination  with  Complete  Pivoting 

To  effect  a  complete  pivotal  strategy  involves  examining  all  the  remaining 
coefficients  for  the  largest  at  each  pivotal  stage  and  exchanging  rows  and 
columns  to  move  that  element  to  the  pivotal  position.  In  our  system  the  lat 
gest  coefficient  is  4,  the  coefficient  of  in  equation  two.  Exchanging 

equation  two  with  equation  one  and  column  three  with  column  one  we  obtain: 

4X3  4  0  t  Xx  =  -2 

2X3  +  X2  4  XT  =  1 

-2X3  4  X2  4  2X1  =  7 

This  achieves  the  desired  result  of  getting  the  largest  coefficient  in  the 
first  pivotal  position,  but  we  have  had  to  abandon  our  matrix  notation,  as 
switching  columns  would  have  lost  the  identities  of  the  dependent  variables. 
As  before,  a)  dividing  4  into  2  and  multiplying  the  first  equation  by  .5  and 
subtracting  the  first  equation  from  the  second;  and  b)  dividing  4  into  -2  and 
multiplying  the  first  equation  by  .5  and  subtracting  the  first  equation  from 
the  third  we  obtain: 


EBJLfJJ  MJ  ww  A  Ul1  Jl  A  VJr/wiU’wviw 


4X3  +  0  + 


Xj  =  -2 


X2  +  . 5Xj  =  2 


X2  +  2.5Xj 


Continuing  our  complete  pivotal  strategy  we  examine  equations  two  and  three 
and  find  we  must  switch  them  and  also  switch  columns  two  and  three  to  place 
2.5,  the  largest  coefficient  in  the  second  pivotal  position  obtaining: 


4X3  ♦ 


Xj  +  0  =  -2 


2.5Xj  +  X2  =  6 


.5X1  +  X2  =  2 


Dividing  2.5  into  .5  and  multiplying  equation  two  by  .2  and  subtracting  the 
result  from  equation  three  yields 


4X3  + 


X1  +  0  =  -2 


2.5Xj  ♦  X2  =  6 


.8X2  =  .8 


By  back  substitution  we  again  obtain  X2  =  l,  X^=2,  and  X3=-l. 

Before  we  justify  this  somewhat  tiresome  detailed  explanation,  note  that 
we  went  about  our  elimination  by  dividing  the  pivot  into  the  element  to  be 
eliminated,  multiplying  the  pivot  row  by  that  quotient  and  subtracting  the 
result  from  the  row  we  were  operating  on.  This  is  important.  We  are  seeking 
the  most  efficient  algorithm  possible--a  numerical  method  that  is  the  best  of 
all  possible  methods,  taking  the  least  amount  of  steps  and  producing  the  most 


accurate  results,  given  the  conditions  of  the  problem.  Notice  we  are  not 
calculating  A  *  in  finding  X.  We  hardly  ever  need  A  What  we  want  is  A  *b. 


Therefore,  in  our  algorithm  we  are  going  to  factor  A  into  two  triangular 
matrices  L  and  U  such  that  LU=A.  L  will  be  lower  triangular  witn  l's  on  the 
diagonal,  zeros  above.  The  elements  below  the  diagonal  will  be  precisely 
those  quotients  obtained  by  dividing  the  element  in  that  position  by  the  pivot 
element  in  that  column.  U  will  be  the  upper  triangular  matrix  we  obtained 
prior  to  starting  the  back-substitution  with  the  pivot  elements  on  the  diago¬ 


nal,  and  zeros  below.  This  approach  not  only  allows  us  to  solve  any  system  of 


n  equations  in  n  unknowns  in  n  /3  operations,  but  the  solution  X'  for  any  new 

2  2 

right  hand  side  b'  can  be  found  in  only  n  operations  ,  allowing  us  to  afford 


extensive  iterative  gambits  with  the  vector  b  or  the  coefficients  A.  The  LU 
decomposition  of  A  in  our  first  elimination  example  is: 


a 


j* 


k 


8 


The  argument  for  the  extremely  conservative  complete  pivotal  strategy  is 
somewhat  tenuous;  one  cannot  assert,  as  in  the  case  of  Gaussian  elimination, 
that  it  is  the  most  efficient  of  all  general  strategies  for  a  high-speed, 

digital  machine.  The  most  obvious  (and  trivial)  benefit  is  that  it  automatic 
ally  avoids  having  to  deal  with  a  zero  pivot.  One  always  has  a  value  in  the 
pivotal  position  if  any  coefficients  greater  than  zero  remain.  This  has  to  be 
dealt  with  in  any  event  and  complete  pivoting  does  it  de  rigueur.  The  product 
of  the  diagonal  elements  of  U  is  the  determinant  of  the  matrix  A.  If  a  pivot 
is  extremely  small  (close  to  0),  the  determinant  is  close  to  0,  possibly  indi¬ 
cating  the  system  is  ill  conditioned,  has  a  singularity,  or  is  unstable. 
These  problems  can  be  dealt  with  individually  or  collectively  by  inspection, 
but  for  an  automatic  algorithm  for  all  seasons,  choosing  the  largest  pivots 

gives  one  the  most  reasonable  chance  for  an  accurate,  stable  solution.  The 
goal  is  to  try  to  keep  the  pivots  as  near  the  same  magnitude  as  possible.  We 
can  do  no  better  than  quote  Professor  Wilkinson,  the  undisputed  authority  on 
the  subject:  "We  do  not  claim  that  selecting  the  largest  elements  as  pivots 

is  the  best  strategy  and  indeed  this  will  often  be  far  from  true.  We  do  claim 

that  no  alternative  practical  procedure  has  yet  been  proposed ...  complete 

3 

pivoting,  on  the  other  hand,  is  never  a  poor  strategy."  A  final  comment  on 
this  admittedly  sketchy  argument  for  complete  pivoting  has  to  do  with  floating 


10 


point  computations.  Uhen  numbers  are  represented  in  machines  by  a  mantissa  m 

multiplied  by  10  raised  to  some  power  x  (n=ml0  where  . 1 < | m | < 1 ) ,  obviously, 

adding  or  subtracting  two  numbers  whose  exponents  differ  say,  by  a  factor  of 

two,  causes  round-off  error  (.246  +  .00193  =>  .247)  if  the  number  of  digits  of 

the  numbers  m  are  of  the  order  of  the  word  length  of  the  machine.  As  there 

are  one  third  of  a  million  operations  in  the  solution  of  a  100  x  100  matrix 
3 

(n  /3  steps),  round-off  error  can  become  a  non-trivial  factor  in  the  solution 
if  extremely  small  pivots  arise  in  the  course  of  the  computation.  It  is, 
therefore,  wise  to  minimize  the  number  of  operations  and  maximize  the  size  of 
the  pivot  elements  if  using  floating  point  numerical  operations. 

We  return  now  to  the  solution  of  the  linear  system  (2)  representing  the 
regression  model  of  our  surgical  tent  problem  which  is  accomplished  by  GAUSS, 
the  program  shown  in  figure  1.  First,  notice  the  services  afforded  the 
analyst  by  FM.  We  DO  ‘ZIS  to  allow  us  to  send  the  result  of  our  efforts  to 
any  addressable  device,  we  set  DUZ(0)="@"  to  meet  security  constraints  on  our 
files,  and  finally  we  DO  “DICRW  to  have  FM  fetch  the  number  of  equations  (NQ) 
and  the  file  number  (FN)  when  we  respond  with  4003  to  the  "OUTPUT  FROM  WHAT 
FILE?"  query.  The  number  of  variables  (NV)  are  then  ascertained  by  addressing 
the  data  dictionary  by  $P('DD(FN,  0),"A",4);  we  subtract  2  from  this  number  as 
we  do  not  count  the  name  of  the  equation  or  the  inhomogeneous  terms  as  varia¬ 
bles  (see  equation  (4)).  Admittedly,  one  must  know  what  information  exists  in 
the  tables  represented  by  the  'DD  and  ‘DIC  globals  but  this  is  exactly  analo¬ 
gous  to  knowing  what  information  we  put  in  our  tables  in  *DIZ,  the  user-crea¬ 
ted  files.  (Incidentally,  all  variables  beginning  with  T  in  GAUSS  are  tempoi 
ary  such  as  TXP;  the  rest  of  the  variables  are,  hopefully,  self-explanatory.) 

We  create  our  working  version  of  Ax=b  by  loading  X(J,K)  witli  the  coeffi¬ 
cient  values  and  the  inhomogeneous  terms  from  ‘DIZ(FN, J ,0) .  X(J,K)  is  out 
relation  N0RMEQ  (File  4003)  with  domains  K,  range  J. 

At  this  point,  we  give  each  instance  in  the  relation  a  unique  name  to 
enable  us  to  switch  columns  without  losing  original  domain  identity--we  give 
each  attribute  occurrence  X(J,K)  the  name  NV*(J-1)+K  which  translates  to  an 
integer  which  we  place  in  XP(J,K).  Consider  a  3x3  matrix  X  with  the  values 
X(J,K);  if  we  placed  these  elements  in  a  string--as  the  computer  does--we 
have,  with  the  elements  numbered: 


GAUSS  ;DRH/NHRC  ;  12/1/87  ;  COMPLETE  PIVOTAL  GAUSSIAN  FORWARD  ELIMINATION 
D  *%ZIS  S  DUZ(0)="@" , U=" * "  D  ' DICRW  S  NQ=%, FN=$P( Y , " ' " ) 

S  NV=$P( "DD(FN ,0),"'",4)-2  I  $D(‘DD(FN, .001,0))  S  NV=NV-1 
F  J=l: 1:NQ  F  K=1:1:NV  S  ( LU( J , K) , U( J ,K) , L( J ,K) )=0 
F  J=1 : 1 :NQ  F  K=1:1:NV+1  S  X( J , K) =$P(* DIZ( FN, J , 0) , " ‘ " , K+l ) , XP( J , K) = 
NV*( J-l )+K 

MAIN  F  PIV=1 : 1 :NQ-1  D  ABSMAX 

F  J=1:1:NQ  F  K=1:1:NV  S  U( J ,K)=X( J ,K) 

F  J=1 : 1 :NQ  F  K=1:1:NV  I  J=K  S  L(J,K)=1 
G  RESULTS 

ABSMAX  S  XMAX=$$MAX(PIV,NQ,NV,.X) 

COLSW  D  CSW(PIV,NQ, .X) 

ROWSW  D  RSW(PIV,NV, .X, .L) 

PIVOTS  F  J=PIV+1 : 1 :NQ  S  L( J, PIV)=X( J , PIV)/X(PIV, PIV) 

F  J=PIV+l:lsNQ  F  K=1:1:NV+1  S  X( J ,K)=X( J , K) -( L( J , PIV)*X( PIV , K) ) 

Q 

RESULTS  F  K=NV :  -1 : 1  S  ANS(K)=X(K,NV+1)/X(K,K) 

F  K=NV-1 : -1 : 1  F  J=NV:-1:2  S  ANS(K)=ANS(K)-(ANS(J)*X(K, J)/X(K,K) ) 

Q: J-K=l 

F  K=NV:-1:1  D  HASH(XP(K,K) ,NV, .ROW, .COL)  S  NAME(K)=$P( ‘ DD( FN , COL, 0) , 

F  K=1 : 1 :NV  W  !!, "VARIABLE :  " , NAME(K) , "  =  " , $J( ANS(K) ,  9 , 6) 

LI  D  MM( .L, .U, -LU)  F  J=1:1:NQ  W  !  F  K=1:1:NV  W  $J(LU(J ,K) ,8,0) 

QUIT 

; SUBROUTINES :  ROWSWITCH , COLSWITCH , MAXELEMENT , MATRIXMULT , HASH 
RSW(P,K,X, L)  NEW  T1,TXP,TRX 

F  Tl=l : 1 : K+l  S  TXP( P , T1 )=XP( P, Tl ) , TRX( P ,T1 )=X( P, T1 ) , X( P ,T1 )= 

X(TJ,T1),XP(P,T1)=XP(TJ,T1),X(TJ , T1 )=TRX(P ,T1),XP(TJ,T1) =TXP( P ,T1 ) 
F  Tl=l : 1:NQ  S  TRL(P,T1 )=L(P ,T1 ) , L(P,T1 )=L(TJ ,T1 ) , L(TJ ,T1 )=TRL(P, Tl ) 

Q 

CSW(P,J,X)  NEW  Tl ,TXP,TCX 

F  Tl-lsls J  S  TXP(T1 ,P)=XP(T1,P) ,TCX(T1 ,P)=X(T1,P),X(T1,P)=X(T1,TK), 
XP(T1,P)=XP(T1,TK),X(T1, TK)=TCX(T1 , P) , XP(T1 , TK) =TXP(T1 , P ) 

Q 

MAX(P , J , K , X)  NEW  MX,T1,T2,TX 
S  MX=-9999999999 

F  T1=P: 1 : J  F  T2=P:1:K  S  TX(T1,T2)=X(T1,T2)  I  TX(T1,T2)<0  S  TX(Tl,T2) 
-TX(T1 , T2 ) 

F  T1=P: 1 : J  F  T2=P:1:K  S:TX(T1 ,T2)>MX  MX=TX(T1 ,T2) ,TJ=T1 ,TK=T2 
QUIT  K*(TJ-1 ) +TK 
MM(X, Y, Z)  NEW  M,N,P,R,I,J,K,V,W 

F  J=1 : 1  Q: ' $D(X( J , 1 ) )  F  K=0:1  S  V=$0(X(J,K))  Q:V="" 

S  M=J-1 ,N=K 

F  J  =  1 : 1  Q:  '  $D(  Y(  J ,  1 ) )  F  K=0:1  S  W=$0(Y(J,K))  Q:W=-"" 

S  R=J-1 , P=K  I  N'=R  W  !,"#  COLS  OF  A  MUST  EQUAL  #  OF  ROWS  OF  B"  Q 
F  1=1 : 1 :M  F  J=1:1:P  F  K=1:1:N  S  Z( I , J ) =Z( I , J )+ ( X( I , K)*Y(K , J ) ) 

Q 

HASH(POS,NV,R,C)  I  POS<NV! (POS=NV)  S  R=1,C=P0S  Q 
I  P0S#NV=0  S  R=POS/NV , C=NV  Q 
E  S  C=P0S#NV,R=(P0S-C)/NV+1  Q 

Figure  1.  Gaussian  Elimination  Algorithm 
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X(l, 1) ,  X(l,2),  XI, 3) ,  X(2, 1) ,  X(2,2.) ,  X(2,3),  X(3,l)t  X(3t2),  X(3,3)  (6) 
Row  1  Row  2  Row  3 

NV  (the  number  of  columns  of  the  matrix)  times  J  (the  row  of  the  element  we 
are  naming)  gives  us  the  maximum  value  the  name  can  take — if  we  subtract  NV, 
the  number  of  elements  in  the  row,  we  obtain  the  name  of  the  last  element  in 
the  previous  row — now  adding  K  gives  us  the  name  (or  position  if  you  prefer) 
that  uniquely  identifies  X(J,K).  The  linear  position  of  X(3,2),  for  example, 
is  3  times  3  (its  maximum  value  possible)  minus  the  number  of  elements  in  the 
row  (3)  yielding  the  position  of  the  last  element  in  the  previous  row,  plus  2 
(its  position  in  the  row  under  consideration)  to  yield  8.  To  undo  this  is 
slightly  more  complex.  The  reader  is  directed  to  the  subroutine  HASH  where 
the  original  column  is  ascertained  by  manipulating  the  number  of  variables 
(NV)  and  POS(ition)  (value  in  XP(J,K)). 

We  now  start  by  finding  the  largest  X(J,K)  (in  an  absolute  sense)  by 
defining  the  extrinsic  variable  $$MAX  and  calling  subroutine  MAX  which  returns 
the  position  name  and  places  TJ  and  TK  in  the  symbol  table  (these  are  the 
temporary  J,K  of  the  largest  element  of  the  matrix  on  the  first  pivotal  pass). 
The  goal  being  to  place  this  element  in  X(PIV,PIV),  ve  now  switch  columns, 
passing  the  entire  X(J,K)  matrix  to  the  subroutine  CSV  by  the  call  by  refer¬ 
ence  .X  .  As  each  column  element  is  manipulated,  its  name  XP(J,K)  gets 
identical  treatment.  The  instruction  lines  in  both  CSV  and  RSV  are  simple  but 
confusing  due  to  the  intermediary  temporary  positions  (TCX,  TRX,  and  TXP) 
necessary  to  etfect  the  bookkeeping.  Notice  in  Row  that  the  elements  in  our 
lower  triangular  matrix  L  are  also  shifted  to  maintain  their  original  func¬ 
tional  integrity  (i.e.  LU=A). 

Once  we  have  MAX  in  the  first  pivotal  position,  we  do  the  forward  elimina¬ 
tion  steps  in  MAIN+7  and  MAIN+8.  The  values  of  L  are  computed,  and  the  cur¬ 
rent  matrix  X(J,K)  has  the  proper  elements  eliminated.  Ve  then  go  back  to 
MAIN,  inclement  PIV  by  1  and  repeat  the  process  until  Vi  have  'J,  the  upper 
triangular  matrix  we  need  for  back-substitution. 

The  back-substitution  begins  at  RESULTS  where  we  solve  for  the  first 
unknown-dividing  the  inhomogeneous  term  X(K,NV-^1)  by  variable  X(K,K):  in  our 
present  case  ANS(3)  =  X(3, 4 )/X( 3, 3) .  It  should  be  clear  at  this  point  we  have 
no  idea  which  variable  it  is,  but  happily  we  have  its  oLiginal  column  identity 
in  XP(3,3'.  Ve  go  on  to  divide  each  of  the  inhomogenous  terms  by  the  approp  - 
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riate  pivot  (the  diagonal  elements  of  U)  so  that  in  line  RESULTS+1  we  can 
solve  for  the  rest  of  the  variables  by  back-substitution. 

We  now  call  our  subroutine  (appropriately  named)  HASH  giving  it  the  origi¬ 
nal  name  of  the  column  of  X(K,K)  which  is  in  XP(K,K)  and  the  number  of  varia¬ 
bles  (NV)  and  ask  it  by  reference  (.ROW,  .COL)  to  place  in  ROW,  COL,  the  J,K 
of  X(K,K)'s  original  position  in  the  matrix  A  of  coefficients.  We  need  only 
COL,  of  course,  but  the  identity  of  the  original  record  could  be  of  interest 
in  other  matters.  Finally,  we  go  to  friendly  FM  and  ask  the  data  dictionary 
to  give  us  the  name  of  the  variable  in  COL  by  SET  NAME(K)=$P( *DP(FN, COL, 0) , 
""",1)  and  print  the  variable  names  and  values  with  line  RESULTS+3.  Line  LI 
of  Gauss  multiplies  the  matrices  L  and  U  to  reproduce  the  original  matrix  A  as 
a  check  of  the  solution.  The  subroutine  MM  first  checks  to  see  if  the  number 
of  columns  of  the  first  matrix  are  equal  to  the  number  of  rows  of  the  second 
matrix  before  attempting  the  multiplication.  This  allows  the  rather  elegant 
matrix  array  passing  notation  MM( . L,  .U,  .LU),  otherwise  one  would  have  to 

pass  the  number  of  columns  of  L  and  the  number  of  rows  of  U,  MM( .  L, .  U, .  LU,C, 
R) ,  which  looks  more  formidable  and  does  not  have  the  face  validity  of  the 
more  succinct  notation.  The  more  responsibility  the  machine  has  in  these 
matters,  the  more  accurate  the  outcome.  The  line  of  code  that  accomplishes 
the  matrix  multiplication 

F  I-ltlxN  F  J=1:1:P  F  K=1:1:N  S  Z(I, J)=Z(I, J)+(X(I,K)*Y(K, J)) 

is  a  reasonable  example  of  MUMPS  programming  parsimony  exploiting  the  power  of 
the  conceptual  model,  the  linearity  of  the  algorithm  and  the  conciseness  of 
the  notation. 

The  solution  proceeded  as  follows,  with  the  original  position  of  the  ele¬ 
ments  shown: 

First  Pivot: 


Original 

Column  Switch 

Row 

Switch 

Elimination 

83,1  25,2  16,3 

25,2 

8,1 

16.3 

87,5 

25,4 

55,6 

87,5 

25,4  55,6 

25,4  87,5  55,6 

87,5 

25,4 

55,6 

25,2 

8,1 

16,3 

0 

.82,1  .20,3 

16,7  55,8  36,9 

55,8 

16,7 

36,9 

55,8 

16,7 

36,9 

0 

.20,7  1.23,9 

Second  Pivot: 

Column  Switch 

Row  Switch 

Elimination 

87,5  55,6  25,4 

87,5 

55,6 

25,4 

87,5 

55,6 

25,4 

0  .20,3  .82,1 

0 

1.23,9 

.20,7 

0 

1.23,9 

.20,7 

0  1.23,9  .20,7 

0 

.20,3 

.82,1 

0 

0 

.79,1 

Back-substituting  yields 

Xl=20191.67,  X2=4133.33,  X3=758.33 


1 

0 

0 

87 

55 

25 

.63218 

1 

0 

U  = 

0 

1.23 

.1954 

.28736 

.15888 

1 

0 

0 

.78505 

87  55  25 

8  25  16 

A  =  LU  = 

55  36  16 

= 

25  87  55 

25  16  8 

16  55  36 

Y'  =  20192.67  +  4133.33  VI  +  758.33  V2 

which  says  that  our  best  estimate  is  that  each  operating  table  adds  $4133  to 
the  cost  of  the  surgical  tent  and  each  extra  light  adds  $758. 

Discussion 

The  basic  ingredients  for  problem  solving  with  a  computer  are  the  compu¬ 
ter,  the  instruction  set  to  run  the  computer,  and  some  data.  Successfully 
entering  information  represented  by  data  into  a  computer  brings  the  analyst 
face  to  face  with  three  fundamental  issues:  (a)  where  are  the  data  and  how  do 
we  get  there;  (b)  what  are  the  data;  and  (c)  what  do  we  do  with  it  (assuming 
we  got  there!  "There  is  nothing  more  embarrassing  to  a  database  administrate: 
than  being  asked  if  his  database  contains  specific  information  and  his  reply 
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ing  after  a  week's  examination  of  the  database  that  he  does  not  know."  ) 

The  computer  machinations  are  themselves  linear--rows  of  is  and  Os  at  one 
conceptual  level.  We  instruct  the  computer  to  sort  these  finite  sequences  and 
we  eventually  compare  the  final  arrangement  to  some  external  criteria  to  judge 
the  value  of  the  computed  information.  The  models  with  which  we  choose  to 
represent  reality  and  those  with  which  we  are  most  comfortable,  are  in  the 
main,  linear.  If  the  data  are  discrete,  as  most  of  the  data  in  our  databases 
are,  we  have  all  the  power  of  linear  algebra  at  our  disposal.  If  we  choose  to 
display  our  data  in  two  dimensional  "flat"  (no  pointers  -each  cell  has  a 
discrete  value)  tables,  we  have  the  power  of  set  theory  to  create  unions, 
intersections,  differences,  projections,  selections,  divisions,  and  joins  of 
our  data. 


What  we  strive  for  in  storing  data  is  a  consistent,  non-redundant  struc¬ 
ture  that  permits  access  to  every  datum.  If  the  structure  is  independent  of 
applications  used  to  derive  the  information  from  the  data,  great  economy  is 
realized  in  creating  applications,  maintaining  the  database,  and  in  accessing 
the  information  by  others.  Of  the  techniques  (hierarchical/networking  and 
relational)  currently  in  fashion  to  store  data,  a  relational  database  manage¬ 
ment  system  has  the  most  promise  because  of  its  simple,  but  powerful,  theore¬ 
tical  underpinning.  It  is  said  that  all  mathematics  can  be  embeade  I  in  set 
theory.  This  is  a  strong  inducement  to  use  relational  technology.  Codd's 

rule  "all  information  in  a  relational  database  is  represented  explicitly  at 
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the  logical  level  in  exactly  one  way--by  values  in  tables"  ,  when  adhered  to, 
offers  the  analyst  a  rational  working  environment.  This  rule  demands  that 
"all  information"--table  and  attribute  information,  data  integrity  con¬ 
straints,  and  security  information,  for  example,  are  available  by  simply 

addressing  the  database  with  the  proper  table  name,  row,  and  column. 

The  Veteran's  Administration  File  Manager  falls  short  of  this  stringent 
criterium,  of  course — as  do  all  current  relational  systems--but  FM  has 

incredible  strengths.  FM  "is  used  to  define,  enter,  and  retrieve  information 
from  a  set  of  computer-stored  files,  each  of  which  is  described  by  a  'data 

dictionary'".^  In  FM  a  file  is  a  relation,  a  record  is  a  set,  and  a  field  is 
in  a  domain. 

One  creates  and  defines  a  file  by  using  option  4 — "Modify  File  Attributes" 
in  the  menu  driven  interface.  The  file  attributes  are  stored  in  a  table 
(MUMPS  global)  named  'DIC.  The  field  attributes  defined  are  stored  in  a  table 
named  ‘DD  (the  "data  dictionary").  The  data  are  entered  using  option  1- 
"Enter  or  Edit  File  Entries"  and  stored  in  a  table  named  'DTZ.  The  question 
"where  are  the  data  and  how  do  we  get  there?"  is  answered  by  responding  with  a 
"?"  to  the  prompt  "OUTPUT  FROM  WHAT  FILE?"  in  the  "Print  File  Entries" 
option--FM  lists  all  the  relations,  giving  file  number,  file  name,  and  number 
of  records.  Knowing  the  file  number  you  can  ascertain  what  is  there  by  using 
option  8  "List  File  Attributes"  and  having  FM  print  the  data  dictionary  infor¬ 
mation  for  each  field  in  the  file.  You  can  get  to  any  stored  datum  now  since 
you  know  the  file  and  the  field.  You  simply  SET  DATUM=$P(  ‘DIZ(FILENUMBER. 
RECORDNUMBER , 0 ) , FIELD  NUMBER) . 

The  question  "what  are  the  data?"  was  answered  when  you  created  the  file. 
Each  field  is  allowed  one  of  seven  data  types--date,  numeric,  a  set,  free 


text,  word  processing,  computed  ("virtual"  field),  and  a  pointer.  In  addition 
to  the  data  type,  range  constraints  and  transaction  restrictions  can  be  set  to 
further  ensure  data  integrity .  "What  do  we  do  with  it?"  is  a  function,  of 
course,  of  what  "it"  is.  Multiple  regression  has  to  do  with  numbers,  so  we  do 
arithmetic.  We  exploit  the  linearity  of  the  machine,  of  the  database  model, 
and  the  notation  by  generalizing  the  problem.  We  let  X  stand  for  the  tables, 
J  for  the  records,  and  K  for  the  fields.  We  can  now  manipulate  X(J,K)  as  a 
relation  (set  theory)  or  as  an  attribute  value  (linear  algebra). 

The  case  for  FM  as  a  relational  database  management  system  is  very  strong. 
There  is  no  doubt  in  the  author's  mind  that  the  relational  "view"  is  going  to 
dominate  the  future  in  data  structuring.  FM  has  all  the  utilities  necessary 
to  do  the  job,  and  it  is  just  a  matter  of  time  before  someone  writes  a  struc¬ 
tured  query  language  for  this  monumental  piece  of  public  domain  software  to 
enhance  an  already  incredibly  useful  tool. 

FM  is  written  in  MUMPS,  which  brings  us  back  to  our  original  question:  Is 
MUMPS  a  reasonable  choice  for  numerical  analysis  in  statistical  work? 

One  of  the  bugbears  in  statistical  analysis  is  the  problem  of  missing  data 
leading  to  underestimation  (of  the  mean,  for  example,  when  N  is  taken  as  the 
number  of  sets  and  5)X  does  not  include  values  for  some  members).  Convention¬ 
ally,  the  delimiter  in  MUMPS  globals  is  the  circumflex  ("'"),  defining  what 
exists  between  's  as  an  attribute  value  in  a  relational  database.  "No  data" 
can  take  the  form  of  null  ("*),  the  empty  string  (*""*),  or  the  empty  set 
(‘O').  The  universe  of  possible  attribute  representations  is  a  combination  of 
one  or  more  of  the  128  definitions  (represented  by  the  integers  from  0  to  127) 
in  the  United  States  American  Standard  Code  for  Information  Interchange 
(ASCII)  character  set. 

Consider  this  line  of  MUMPS  code: 

F  J=1:1:NV  S  N(J)=NR  F  K=1:1:NR  I  X( J ,K) ' J$C(0)  S  N(J)=N(J)  1 

where 

J  =  the  number  of  the  domain  (field  number  or  variable,  however  you  are 
picturing  your  data) 

NV  =  the  last  domain  to  be  examined  (the  number  of  variables) 

K  =  the  set  numbei  (record) 

NR  =  the  range  of  J  (the  number  of  records) 

I  =  IF 

X  =  the  relation  (the  table  name  or  file  name) 


]  =  follows,  in  the  sense  of  the  collating  order  of  the  ASCII  character  set 

$  =  signals  what  follows  is  a  special  instruction--evaluate  the 

parenthetical  argument 

C  =  the  symbol  for  the  function  that  takes  the  integer  value  argument  and 

returns  the  ASCII  definition  it  represents  (some  of  which  are  not  printa¬ 
ble). 

0  =  the  decimal  ASCII  value  for  NULL,  the  first  in  the  code  which  ranges  from 

0  to  127 

S  =  "Set"  in  MUMPS 

N  =  the  number  of  sets  to  be  examined. 

This  is  an  existence  theorem.  It  states  "set  N  for  the  field  you  are 
examining  equal  to  the  number  of  records  you  are  going  to  examine.  From  the 
first  record  to  the  last,  look  at  each  field.  If  there  is  no  value  for  that 
record  and  field,  reduce  N  by  one." 

I  X(J,K)' ]$C(0)  says,  if  the  contents  of  X(J,K)  does  not  follow  the  first 

value  in  the  ASCII  code  (which  is  null),  then  the  contents  must  be  null--if  it 

ain't  in  the  universe,  it  doesn't  exist!  (Remember,  by  definition,  ""  and  0  do 

exist,  of  course,  and  one  can  evaluate  them  in  whatever  context  they  occur. 

The  empty  string  might  mean  the  subject  refused  comment;  zero  could  be  his/her 

response  to  number  of  children.)  Typically  the  time  to  use  this  powerful 

screen  is  when  one  is  loading  the  working  data  subset  into  resident  memory. 

At  this  time  one  should  also  transpose  the  data  (turn  the  table  on  its 
T 

side--X  (J,K)  =>  X( K , J ) )  so  that  the  rows  consist  of  the  variables  and  the 
columns  are  the  records  in  which  the  variables  occur.  The  discussion  of  the 
existence  rule  followed  the  transposed  concept.  This  is  also  consistent  with 
mathematical  notation  for  a  function  which  is  a  relation  on  coupled  pairs: 

f(x)  =  R(x,f(x))  =  R(x,y)  =  y 

where  x  is  the  domain,  y  is  the  range,  R  is  the  relation,  and  f  signifies  a 
function  or  rule  on  x  (f  s  R).  It  takes  three  times  as  long  to  sum  the  ele¬ 
ments  of  a  matrix  (X(J,K))  by 

F  K=l:l:NCols  F  J=l:l:NRovs  S  S=S+XT(J,K) 
than  it  takes  doing 

F  K=1 : 1 :NCols  F  J=l:l:NRovs  S  S=S+XT(K,J) 

on  a  VAX  11/785.  Essentially,  what  the  second  line  amounts  to  is  reading  all 
the  descendants  of  a  node  before  preceding  to  the  next  node.  If  one  does  not 


transpose  the  data,  the  machine  has  to  repeatedly  skip  from  node  to  node  for 
each  datum  (see  the  linear  array  in  equation  (6)).  An  analogy  is  the  process 
of  shopping  in  a  supermarket.  You  go  down  the  aisles  (records)  and  pick  from 
the  bins  (fields).  You  seldom  shop  aisle  1  bin  1,  then  aisle  2  bin  1  to  aisle 
12  bin  1  and  then  start  over  at  aisle  1  bin  2,  etc.  Unfortunately,  poor  pro¬ 
gramming  strategies  like  this  occur  all  too  often  in  the  misuse  of  MUMPS  and 
lead  to  opinions  MUMPS  is  slow  and  unsuitable  for  scientific  work.  The  penal¬ 
ties  are  severe  in  local  memory  but  if  extended  to  reading  globals  the  in¬ 
crease  in  disk  reading  time  is  astronomical.  (The  ratio  of  the  two  approaches 
remains  3:1,  but  the  time  goes  from  milliseconds  to  minutes.) 

The  power  and  conciseness  of  the  MUMPS  language  is  without  parallel  in  the 
routine  for  deriving  the  square  root  of  positive  real  numbers: 

S  Y=0  Q:X'>0  S  Y=l+X/2 
L  S  T=Y, Y=X/T+T/2  G  L: Y<T 
and  in  this  Shell  sort  (ascending  order)  routine^: 

S  G=N\2 

S2  F  I=G+1 : 1:N  S  L=I-G  I  (Y(J, L)>Y(J,L+G))  S  T=Y(J,L),Y(J,L)=Y(J, 

L+G) , Y( J, L+G)=T  G  S2 
S  G=G\2  Q:G' >0  G  S2 

which  divides  a  list  Y  of  N  items  into  two  halves,  compares,  in  order,  the 
elements  of  the  first  half  with  the  elements  of  the  second  half,  and  switches 
them  if  the  first  is  greater  than  the  second.  After  the  first  ordering  the 
list  is  divided  into  quarters  and  the  process  repeated;  finally  adjacent 
elements  are  compared  completing  the  ordering. 

These  examples  strongly  support  the  contention  MUMPS  can  be  used  procedur- 
ally  in  a  numerical  context.  If  the  analyst  clearly  differentiates  between 
storing  things  and  manipulating  them,  efficient  algorithms  can  be  written  in 
MUMPS.  File  structuring,  sensible  blocking  on  disks,  emphasis  on  exploiting 
the  size  of  the  partition,  and  structured  programming  will  do  a  lot  to  dispell 
the  notion  MUMPS  is  slow  when  manipulating  numbers. 

The  emphasis  on  networking  in  the  sixties  stemmed  from  the  electrical 
engineers  who  were  our  first  programmers.  This  influence  is  seen  in  the  idea 
of  embedding  data  in  the  subscript  of  a  MUMPS  global.  The  path  dependencies 
resulting  have  proven  to  be  almost  intractable  for  analyst;.  The  current 
progress  in  compiling  MUMPS  lends  encouragement  to  those  writing  code  that  is 
linear  in  nature,  thus  producing  efficient  object  code. 


\\ 


The  fact  that  ISM  utilizes  19  fixed  point  decimal  digits  on  DEC  machines 
is  a  numerical  analyst's  dream — this  is  automatic  double  precision — when 
compared  to  the  round-off  error  nightmares  the  fifties  and  early  sixties  with 
eight  and  ten  digit  floating  point  machines.  String  manipulation  features 
like  $F  and  $L  allow  incredible  shortcuts  to  examining  numbers,  cutting  the 
code  necessary  to  write  a  logarithm  routine,  for  instance,  to  a  dozen 
lines  —  that  are  as  fast  as  all  but  the  most  strongly  typed  languages  like 
PASCAL — certainly  comparible  to  FORTRAN  if  the  I/O  overhead  is  considered. 
(This  is  intuitive;  benchmarking  comparisons  would  be  interesting,  using 
compiled  MUMPS.) 

The  future  of  MUMPS--its  longevi ty--is  assured.  The  marriage  of  PROLOG 
with  MUMPS  in  GNOSIS^  is  based  on  MUMPS  string  manipulation  capabilities, 
global  aata  structure,  and  indirection  facility.  The  use  of  MUMPS  in  database 
management  is  a  function  of  its  subscripting  and  efficient  global  structuring. 
The  parameter  passing  facilities  now  available  make  possible  maximally  compact 
code.  MUMPS  is  a  language  of  the  future  if  its  strengths  are  exploited  and 
not  denigrated  by  unintelligent  use. 
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Figure  2.  Bi-Variate  Normalized  Gaussian  Distribution 

Conclusion 

The  simplicity  of  the  relational  data  structure  can  lead  to  powerful  ana¬ 
lytical  processes.  it  one  can  visualize  figure  2  as  a  K-J  plane  with  the 
value  of  the  X(K,J)  attribute  represented  by  the  height  of  the  curve  above  the 


point  K,J,  one  has  a  three-dimensional  view  of  the  data  distribution.  Figure 
2  happens  to  be  the  bi-variate  normal  Gaussian  distribution,  but  distribution 
functions  can  be  fitted  to  real  data  in  columns  of  relational  tables  and 
easily  graphed  yielding  dramatic  insight  into  the  data  relations  and  distribu¬ 
tions  over  the  range  of  the  domains. 

Two  roles  for  the  symbols  K  and  J  are  being  utilized:  (1)  K  and  J  as  a 
coupled  pair  (of  integers,  for  example)  to  locate  a  value  of  a  point  in  a 
plane;  (2)  K  and  J  as  the  place  holders  for  the  name  of  a  column  in  a  table 
and  the  name  of  a  row  in  that  table. 

If  we  designate  the  table  as  X,  we  can  write  X(K,J)  to  describe  the  table 
X  with  domain  K,  range  J.  This  now  allows  us  to  treat  the  X  in  set  theoretic 
fashion  employing  all  the  operations  of  that  branch  of  mathemat ics--i t  does 
not  matter  what  the  K  is  at  this  level  as  long  as  some  K  in  the  set  we  are 
operating  on  are  like  some  K  in  other  sets  (for  certain  operations). 

If  we  assign  cardinal  numbers  to  K  and  J,  we  can  manipulate  the  values 
that  each  pair  represent  algebraically,  as  we  did  in  GAUSS,  arriving  at  num¬ 
bers  called  statistics  which  purport  to  say  something  about  how  the  attributes 
in  a  domain  are  related  to  some  cri ter ia~-say  the  mean  of  all  the  salaries,  or 
predict  the  contribution  of  different  elements  to  the  price  of  a  surgical  tent 
as  a  linear  combination  of  the  domains  involved  when  those  domains  are  com¬ 
bined  to  form  a  basis  for  a  vector  space  (least  squares  solution  to  an  over¬ 
determined  system  of  equations). 

The  relational  database  structure  is  sound  mathematically.  It  allows  us 
to  gracefully  handle  issues  of  data  integrity,  independence,  and  access.  When 
viewed  as  presented  here,  it  can  be  exploited  with  all  the  power  of  set  theory 
and  linear  algebra,  which  in  turn  leads  to  fast,  efficient  machine  code. 

Theoretical  Implications 

There  are  at  least  three  emerging  areas  that  the  use  of  MUMPS,  relational 
technology,  and  the  underlying  mathematical  calculus,  are  eminently  suited. 
The  first,  in  conceptual  simplicity,  is  distributed  parallel  processing.  It 
is  not  difficult,  analytically  or  from  the  hardware  standpoint,  to  visualize 
two  or  more  80386  processors  operating  on  the  columns  of  matrices.  One  could 
compute  every  conceivable  statistic,  probability  density  function,  or  sort 
almost  instantaneously  on  each  column. 


The  second  theoretical  consideration  is  the  use  of  MUMPS  and  set  theory  in 
knowledge  based  systems — so  called  artificial  intelligence.  The  indirection 
feature  of  MUMPS  that  allows  the  program  to  write  code  modifying  itself  is  of 
enormous  import  in  problems  involving  servo  mechanisms,  cybernetics,  analysis, 
and  decision  making  (there  is  no  free  lunch;  this  same  facility  exacts  a  price 
in  the  compilation  of  the  code).  The  string  manipulating  features  of  MUMPS 
easily  duplicate  LISP  features  with  appropriate  coding.  The  local/global 
storage  structure  of  MUMPS  reduces  I/O  to  child's  play  and  is  maximally  effi¬ 
cient  in  arrangement  and  space.  The  choice  of  MUMPS  by  the  Japanese  to  drive 
PROLOG  in  GNOSIS  in  creating  a  fifth  generation  language  was  not  whimsical. 

The  third  body  of  ideas  beginning  to  create  widespread  interest  is  the 

Q 

paradigm  of  computational  reflection  .  This  concept  is  illustrated  in  this 
paper  in  the  programming  to  keep  track  of  the  original  domain  of  the  elements 
of  the  matrix  A.  Another  example  is 

P  J=l:l  Q:'$D(X(J,1))  P  R=0:1  S  V=$0(X(J,K))  Q:V=""  S  M=J-1,N=K 
in  the  subroutine  MM  where  we  ask  the  machine  to  determine  if  X(J,K)  is  a 
proper  object  (i.e.  checking  to  see  if  the  number  of  columns  of  matrix  A  equal 
to  the  number  of  rows  of  matrix  B)  for  the  operation  to  follow  and  to  inform 
us  if  it  is  not.  This  might  be  thought  of  as  meta-programming.  It  is  a 
program  about  the  program.  These  examples  illustrate  exactly  what  an  object 
oriented  programming  schema  would  accomplish — we  have  linked  meta-information 
about  the  object  to  the  manipulation  of  the  object.  To  be  able  to  write 
commands  based  on  what  has  happened  in  the  meta  phenomena,  as  we  can  in  MUMPS 
by  using  indirection,  allows  the  programmer  an  unprecedented  level  of  control. 

Conceptually,  MUMPS  globals  are  inherently  relational.  The  Veteran's 
Administration  File  Manager  is  essentially  a  collection  of  information  in 
MUMPS  globals.  Relational  database  structures  are  two  dimensional  tables  that 
are  most  efficiently  manipulated  by  set  theory.  The  entities  of  the  tables 
are  amenable  to  the  rigors  of  linear  algebra.  The  resolution  of  a  classical 
problem  in  numerical  linear  algebra  demonstrates  the  ease,  simplicity,  and 
efficiency  of  the  use  of  MUMPS,  FM,  and  relational  techniques. 

The  choice  of  a  computational  language  that  contains  the  tools  needed  for 
innovative  research  in  a  mathematically  sound  manner  is  prudent.  MUMPS,  FM, 
and  relational  concepts  offer  that  virtue. 


Author's  Note:  Professor  Strang's  delightful  book  "Linear  Algebra  and  Its 
Applications"  was  relied  upon  heavily  in  this  paper.  The  numerical  example 
for  Gaussian  elimination  is  from  that  book  and  was  used  for  ease  of  verifica¬ 
tion.  Professor  Wilkinson's  two  classics  "The  Algebraic  Eigenvalue  Problem" 
and  "Rounding  Errors  in  Algebraic  Processes"  served  as  authority  to  support 
the  remarks  on  the  process  and  to  provide  examples  to  check  the  error  bounds 
on  the  algorithm.  The  author  lays  modest  claim  to  the  hashing  stratagem  as  an 
original  contribution. 
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